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Q ■ Abstract 

I Fermionic hnear optics is a hmited form of quantum computation which is known to be 

efficiently simulable on a classical computer. We revisit and extend this result by enlarging 
Q^. the set of available computational gates: in addition to unitaries and measurements, we 

, K I allow dissipative evolution governed by a Markovian master equation with linear Lindblad 

Ch ' operators. We show that this more general form of fermionic computation is also simu- 

^ ■ lable efficiently by classical means. Given a system of fermionic modes, our algorithm 

Q^i simulates any such gate in time 0{N^) while a single-mode measurement is simulated in 

time 0{N'^). The steady state of the Lindblad equation can be computed in time 0{N^). 

> ; 

^ ; 1 Introduction 

^ \ Quantum dynamical processes have inspired numerous models of computation: adiabatic quan- 
C^l I tum computation [H], dissipative quantum computation |29l |19] or computations based on 
^ ' modular functors [11] are all examples of computational models associated with certain time- 
^ ■ evolving quantum systems. In trying to characterize their computational power, arguably the 

>■ . most practically relevant question is how they compare to universal classical respectively quan- 
^ \ tum computation. Indeed, if computations in a model are efficiently simulable on a classical 
\ computer, the corresponding physical system may be accessible to numerical studies, but is 
unlikely to be a suitable substrate for building a quantum computer. In contrast, simulability 
by quantum circuits means that the underlying physics could be studied using a quantum com- 
puter, the prime application of such machines originally envisioned by Feynman [9J. Finally, in 
cases where the computational resources provided by a model are sufficient to implement univer- 
sal quantum computation, the corresponding physical system is a candidate for the realization 
of a quantum computer. 

A number of physically motivated models can be understood as the result of restricting the 
available set of initial states, gates and measurements in the standard quantum circuit model. 
For topological [HI [221 110] or permutational [13] quantum computing, there is a preferred 
initial (vacuum) state, the available gates represent braid group generators or transpositions 
and there is a set of allowed (charge) measurements. In bosonic linear quantum optics [2T] , the 
available repertoire includes preparation of the vacuum initial state, single photon sources, beam 
splitters, phase shifters and photo-detectors. In fermionic quantum optics [261 120] , we permit 
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preparation of the vacuum state, free unitary evolution and occupation number measurements. 
The computational power of topological computing depends on the representation: for example, 
it is universal in the case of the Fibonacci model but classically simulable for Ising anyons. 
Bosonic linear optics was shown to be universal for quantum computation. Fermionic 
linear optics can be simulated efficiently on a classical computer [261 120] (see below). 

The model of dissipative quantum computing is of a conceptually different origin: it is the 
result of considering more general quantum dynamics beyond unitary evolutions. Here the 
standard unitary gate set is augmented or completely replaced by 'dissipative gates'. Mathe- 
matically, the latter are completely positive trace-preserving maps corresponding to the time- 
evolution under a Markovian master equation. Such Markovian evolution is usually associated 
with noise in implementing a quantum computer, and one seeks to reduce its detrimental ef- 
fect by the use of error-correction. Recently, however, it was realized that purely dissipative 
processes can actually be useful for quantum computation (see e.g., [121 ESI ED])- In [2S], it 
was shown that dissipative quantum computation is universal for quantum computation. Con- 
versely, Markovian dynamics with local Lindblad operators can be simulated efficiently (i.e., 
with polynomial overhead) on a quantum computer ^19]. 

Here we consider a model of computation that extends fermionic linear optics with dissipa- 
tive processes. We will show that this extended model is still efficiently simulable classically. 
Since our simulation algorithm does not rely on Trotter- type expansions (as e.g., [IS]), dissipa- 
tive processes can be simulated exactly for any evolution time without affecting complexity or 
accuracy. The algorithm reproduces the statistics of measurement outcomes and also provides 
a complete description of the state at any instant in the computation. Additionally, the Liou- 
villian may be non-local, but we restrict to Lindblad operators which are linear in fermionic 
annihilation and creation operators. 

The physics underlying fermionic linear optics - non-interacting fermions - encompasses 
a number of systems of interest in condensed matter physics, including Kitaev's Majorana 
chain [16] or honeycomb model [17]. Such systems exhibit topological order and could be 
used as fault-tolerant quantum memories or topological quantum computers. The classical 
simulation algorithm discussed here may be particularly suited to assess their potential to serve 
such information-processing purposes. Typically, this involves studying not only the dynamics 
of the system, but also the effect of simple manipulations (such as syndrome measurements or 
error correction). As an example, the known simulation technique in the non-dissipative case 
has been used to study the beneficial effect of disorder on the performance of the Majorana 
chain as a quantum memory (5]. The extended toolkit provided here may be applied to evaluate 
the robustness of certain proposals, e.g., for state transfer [321 E3], in the presence of dissipation. 

2 Fermionic linear optics 

Fermionic linear optics is defined in terms of creation and annihilation operators Oj and aj 
satisfying Fermi-Dirac canonical commutation relations {aj,ak} = and {aj,al} = Sj,kl- We 
will refer to such a family of operators {aj, ci]}jLi as A^ Dirac fermions or modes to distinguish 
them from Majorana fermions introduced below (the latter will be denoted by the letter c 
instead of a). The corresponding Hilbert space is spanned by the number states 




where rij G {0, 1} , 



(1) 
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and |0) is the fermionic vacuum state satisfying aj\0) = for all j. A state of the form ([T]) is 
an eigenstate of the occupation number operator Nj = a}-aj with eigenvalue rij. 
The allowed operations defining fermionic linear optics are: 

(i) preparation of the fermionic vacuum |0) 

(ii) measurement of the occupation numbers A'^-, j E S for any subset S of modes, 

(iii) evolution under a quadratic fermion Hamiltonian H for a time t, that is, according to the 
equation of motion 

ip-m,A- (2) 

Fermionic parity preservation implies that if is a linear combination of terms of the 
form Cjajaj giving energy e M to mode j, 'hopping terms' of the form tj,fcajafc + t*j^k'^\(ij 
where tj^^ £ C for j k and 'pair creation/annihilation terms' of the form Sj^fcOj^I + 

These operations can be performed in an arbitrary order, and, in particular, may depend on 
measurement results obtained in the course of the computation. The computation concludes 
with a final measurement whose outcome is the classical result produced by the computation. 

Simulating such a fermionic linear optics computation on a classical computer amounts to 
sampling from the distribution of measurement outcomes at the end of the computation. A 
classical polynomial-time (in N) algorithm for doing so was found by Terhal and DiVincenzo 
and Knill [20] . A stronger form of simulation outputs a description of the state at any instant 
during the computation. An efficient algorithm for this problem was provided by Bravyi [1]. 
We review these techniques in Section [51 

Fermionic linear optics was motivated by the quantum universality of bosonic linear op- 
tics [2T], but is closely related to another computational model: the unitaries ([m]) are so-called 
matchgates. The corresponding computational model - matchgate computation - was earlier 
shown to be classically simulable by Valiant |27j. Josza and Miyake [15] extended the sim- 
ulation results by showing that more general initial (product) states can be allowed without 
losing classical simulability. They also showed that a slight modification of Valiant's gate set 
provides quantum universality (see also [2| for a physically motivated universal gate set extend- 
ing fermionic linear optics). More recently, new and elegant characterizations of the power of 
fermionic linear optics have been provided: Jozsa and coauthors [13] showed that the model is 
equivalent to space-bounded quantum computation (see also [5]), and van den Nest [28] gave a 
characterization in terms of linear threshold gates. 

3 Dissipative fermionic linear optics 

Here we generalize these results as follows: in addition to the unitary set of gates defined 
by dm]), we allow dissipative processes. More specifically, we replace dm]) with 

(iii') evolution under the Lindblad master equation 

j^p = C{p) := -t[H, p] + Y, i^L.pLl - {LlL„ p}) , (3) 
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for a time t, where if is a quadratic fermion Hamiltonian and each Lindblad operator 
is a Unear combination = J2j '^mj'^] + /^M.i'^i with a^j, (3^j G C. 

We show that a fermionic computation composed of (El), (Ell) and (En]) can be efficiently simulated 
on a classical computer. 

Outline 

The remainder of this paper is structured as follows. In Section HJ we discuss fermionic Hamil- 
tonians and Liouvillians and we show that dissipative dynamics with linear Lindblad operators 
preserves the set of Gaussian fermionic states. In Section [5l we review the classical simulation 
of fermionic linear optics. In Section [6l we explain how to simulate dissipative fermionic linear 
optics. 



4 Background 

In this section, we give some background on quadratic open fermion systems with linear Lind- 
blad operators. It will be convenient to work with the Hermitian operators 

C2j-i = ttj + Oj and C2j = i{aj — a]) for j = 1, . . . , . 

The 2N operators ci, . . . , C2n satisfy the canonical anticommutation relations 

{cj, Cfc} = 25j^kl for all j,k = 1, . . . , 2N , 

and will be referred to as Majorana fermions or modes. The Hamiltonian and Lindblad opera- 
tors ([3]) then take the form 

.27V . 2N 

j,k=l j=l 

where c = (ci, . . . ,C2n)- Here H = — G so(2A^, M) is an antisymmetric matrix with real 
entries. In contrast, the vectors £^ = • • • ,C-^,2n) ^ are generally complex- valued. 

4.1 Third quantization description of Liouvillians 

We will present two methods for simulating dissipative dynamics of the form ([3]). One of these 
methods is based on the language of 'third quantization'. This refers to the fact that the 
Liouvillian C is quadratic in a set of fermionic superoperators as discussed by Prosen [23] (see 
also |25j and Dzhioev and Kosov [7j). Here we follow his presentation. 

Let /C be the 4^-dimensional vector space spanned by the monomials c- = c"^ ■ ■ ■ 
where a = (ai, . . . , a2N) G {0, 1}^^. To emphasize the vector space structure, we will write an 
operator p = J2a ^aC- (with G M) as |p)) = J2a '^a\c-} when it is considered as an element 
of /C. The space /C is equipped with the Hilbert-Schmidt inner product ((p|cr)) = 2~^TrpV, 
and with respect to the latter, the monomials |c-)) are an orthonormal basis. The space /C 
decomposes into a direct sum JC = /C~*"©/C~ of spaces spanned by monomials |c-)) with even (/C+) 
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and odd (/C~) parity X]j=i mod 2. In the following, we will often restrict our attention to the 
subspace /C+, although identical arguments apply to /C^. By fermionic parity superselection, 
this is sufficient to cover physical states. 

We define 2N pairs (Sj, aj) of mutually adjoint linear operators acting on /C by 

aj|c-)) = '^a„i|c,0) and a]|0) = 5,^,o|c,0) for j = 1, . . . , 2iV. (5) 

These satisfy canonical Dirac anticommutation relations 

{aj,afc} = {o-j^al) = 5j^k for j, k = 1, . . . ,2N , 

which has motivated the expression 'third quantization'. Observe that the state |c-)) = |J)) cor- 
responding to the identity operator / is the vacuum state associated with this set of operators, 
i.e., aj\I} = for all j. 

It is again often more convenient to work with Hermitian Majorana fermions 

C2j_i = flj + a] and C2j = i{a.j — a}-) for j = 1, . . . , 2N . (6) 

The restriction = £|ac+ of the Liouvillian C to operators supported on /C+ will be called the 
even part of C. Remarkably, this superoperator can be expressed as a quadratic form of the 
operators ([5]) (or equivalently ([6])). More precisely, we have 

= i2+|p)) for all p with supp{p) C JC^ , where C+ = -c-hc , (7) 

and where c = (ci, . . . , C4n) for an antisymmetric complex- valued matrix L = — G so(4A^, C). 
The matrix L depends linearly on H and quadratically on the entries of the vectors specifying 
the Lindblad operators. Explicitly, it is given by [23] 

L2j-i,2fc-i = Hj^fc — + 2Mfcj L2i-i,2fc = 42Mfcj . . 

L2i,2fc = H,- fc + 2M,- - 2Mkj Ls,- 2fc-i = -4iM,- ^ ^ 

with Mj fc = [231 125], it is assumed that the Liouvillian is generic in the sense 

that L is diagonalizable. Here we do not require such an assumption. Note also that in [23], a 
normal form for such matrices L is derived. 

Observe that for any integer n and any state p supported on /C+, we have 

|£!^(p))) = |£+o...o£+(p))) = ^"|p)), 

V ' 

n 

which implies that the superoperator exp(t£+) corresponding to time evolution for a time t 
under is given by 

exp(t£+) = exp(t£+) (9) 

when restricted to operators supported on /C+. Because detexp(A) = exp(Tr(A)), the opera- 
tor Qj is invertible with inverse exp(— 

Finally, consider the adjoint Liouvillian defined by 

£t(0) = z[H, 0] + J2 iKOL, - {LlL,, O}) . (10) 
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Evolution for a time t under C"^ generates the adjoint superoperator of exp(t£), that is, 
Tr(exp(t£) (p)0) = Tr(pexp(t£'l")(0)) for all operators p, O. This follows by comparing the 
derivative 

|Tr(Oe*^(p)) = Tr(0£(p)) = Tr(£t(0)p) = | Tr(e*^^(0)p) 

and observing that the two expressions agree for t = 0. Since exp(tC) is a completely positive 
trace-preserving map, the map exp(t£^) is unital, that is, exp(t£"'')(/) = /. The even part C\_ 
of the adjoint Liouvillian takes the form 

€l = \c-Vc (11) 

where the matrix L^^ = — L* G so(4A^, C) is the Hermitian conjugate of L. 

4.2 Gaussian states 

Here we collect a few facts about Gaussian states of 2N Majorana fermions ci, . . . , C2n (see |1] 
for details and proofs). A Gaussian state p is completely determined by its covariance matrix 

M,-fc = ^Tr(p[c,-,Cfc]) for j,fc = l,...,2iV (12) 

and Wick's formula (see e.g., [U Eq. 17]) 

Tr (pc,,c,, ■ ■ ■ c,J = Pf (M [ji, . . . , j2p]) for all 1 < j i < j 2 < ■ ■ ■ < J2p < 2N . (13) 

Here M[ji, . . . ,j2p] is the submatrix of M of size 2p x 2p obtained by removing all columns 
and rows except those indexed by ji, . . . ,j2p, and Pf denotes the Pfaffian. The antisymmetric 
matrix M G so(2iV, M) can be brought into block- diagonal form by a special orthogonal matrix 

M = i?0 f ° Re S0{2N), Aj G R . 

j=i V J / 

This is called the Williamson normal form of M. The transformation c-^d = Rc, that is, 

2N 



'^ = J2Ra,bCb fora = l,...,2iV 



b=l 



corresponds to the adjoint action of a unitary as explained below. In particular, c[, . . . , Cgjy 
again satisfy canonical commutation relations. We call these the eigenmodes of p. Expressed 
in terms of these operators, p takes the simple form 



P = 7^Ili^ + '^A~A) ■ (14) 
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4.3 Time evolution of Gaussian states 



Evolution of Gaussian states under the dissipative dynamics described by Eq. ([3]) is particularly 
simple due to the following fact. 

Lemma 1. Let C be a Liouvillian for N fermions whose unitary part is given by a quadratic 
Hamiltonian and whose Lindblad operators are linear in the creation- and annihilation operators 
(cf. ([3]) ). Let p be a Gaussian state. Then the time-evolved state e^^{p) is Gaussian for allt > 0. 

Proof. We will construct a family of trace preserving completely positive maps depending 
smoothly on a parameter e > such that each map $e preserves the set of Gaussian states and 



= C(p) (15) 

e=0 



for any state p. This implies that 



p{t) = hm <|.V(p) 



t/r. 



is a limiting point of a sequence of Gaussian states. Since the set of Gaussian states is compact, 
p{t) must be a Gaussian state itself. 

Let us first consider a Liouvillian with a single Lindblad operator, 

£(p) = 2LpLt-{LtL,p}, 

where L is a linear combination of the Majorana operators ci, . . . , C2n with complex coefficients. 
Decompose L = K + iM, where K, M are real linear combinations of the Majorana operators. 
In particular, both K and M are Hermitian. Introduce one ancillary pair of Majorana operators 
C2Ar+i = bi and C2N+2 = &2 representing an environment and consider a unitary operator 

Ue = exp 

Simple algebra shows that 



2e{Kbi + M62) 



U.riU^^ = ri- V2^[Kbi + Mfes, v] - 2eiKbi + Mb2)r]{Kbi + M62) 

+t{{Kb^ + Mb2)\ri} + 0(e=^/2) (^g) 

for any state rj. Define a map 

^,{p) = 1:ieU,PPeUI . (17) 

Here Tr^; represents the partial trace over the environment, and pc is the initial state of the 
environment which we choose as the vacuum state, that is, 

PE = ]^{I - ibib2) . 

Note that the terms in Eq. (fT6!) that contain half- integer powers of e do not contribute to $e(p) 
since Tr^;^! = Tr£;&2 = 0. Using the identities 

TiEbib2PE = i and TiEbipEb2 = -i 
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it is easy to check that 



<l>,(p)=p + e£(p) + 0(e2) 



(18) 



Given a general LiouviUian with muhiple Lindblad operators and the unitary evolution term 
as in Eq. ([3]), the desired map •I'e can be constructed by taking a composition 

$^(p) = e-^-f^^$W o . . . o $('")(p)e^^S 

where ^i^^ is the map defined above with L = L^. By construction, $e is a finite composition 
of unitary evolutions under quadratic Hamiltonians, addition of ancillary two-mode vacuum 
states pe, and partial traces over some pairs of modes. It is well-known that all these operations 
preserve the set of Gaussian states, see e.g. [Ij. Finally, Eq. (fT5|) follows from Eq. (fT8|) and the 
product rule for derivatives. □ 



5 Classical simulation of fermionic linear optics 

In this section, we review the known simulation techniques for a computation composed of the 
operations introduced in Section [2l The basis of the simulation algorithm is the fact 

that the vacuum state |0) prepared by (i) is a Gaussian state and all subsequent operations 
preserve the Gaussian nature of the state. 

The covariance matrix M(0) of the vacuum state |0) at the beginning of the computation 
is given by the non-zero entries 

M(0)2j-i,2i = 1 for J = l,...,2iV (19) 

above the main diagonal. The remaining task is to find update rules for the covariance matrix. 
These updates can be done efficiently as discussed below: measuring (time) complexity in 
terms of the number of additions, multiplications, and divisions on complex numbers that are 
required, measurements and unitary gates can be simulated in time 0{N^). 



5.1 Simulating measurements 

Here we describe the method from [6] which is a more efficient version of Terhal and DiVin- 
cenzo's algorithm (the latter has time complexity O(iV^) when measuring N modes). 

Consider a (non-destructive) measurement of the occupation number Nj = a^aj = ^{I ~ 
ic2j~iC2j) in a Gaussian state p with covariance matrix M. By definition (1121) . the probability 
of obtaining outcome 1 when measuring p is 

P,(l) = Tr(iV,p) = i(l - M2,_i,2,) . (20) 

Let Uj^rij) = (a]aj)"-' (aj-aj)^""^ be the projection corresponding to the measurement outcome 
rij G {0, 1}. The post-measurement state 



p(n,) = "^-^"^^^^^^-^ (21) 
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is Gaussian, as shown in [4J . Its covariance matrix {rij) can be computed from fl20|) and fl2T]) 
using Wick's theorem ( 1T3|1 . giving (cf. [6, Eq. (7.3)]) 

M^^\nj),,, = M,,, - l^J-^M2j-i,,M2,,p + jp^M2j-i,,M,j,p . (22) 

The following algorithm then simulates a measurement of Nj: first, compute (120|) and sample 
a bit nj G {0,1} according to the probability distribution Pr[nj = 1] = -P,(l). Then update 
the covariance matrix according to (122]) . This requires 0{N^) computational steps and one 
(non-uniform) bit of randomness. 

A measurement of a subset S = {ji, . . . ,j\s\} of modes can be simulated by iterative use of 
this procedure, using the recursion relation 

Tr(nj,(njjp(nj, ■■■nj,_J) 

for the probability and post-measurement state after the measurement of the jVth mode, given 
that the measurements of Nj^, . . . , Nj^_^ resulted in the sequence n^^, . . . , nj^_^ G {0, 1}. This 
can be done in time 0(|5| ■ N'^) using \S\ random bits. 



5.2 Simulating unitary evolution 

The classical simulation of unitary dynamics (see PU| [251 [T5]) is particularly instructive for our 
generalization to the dissipative case. Consider time evolution under a quadratic Hamiltonian 
H = -f^(H) (cf. (jlj)) for some time t, starting from a Gaussian initial state p with covariance 
matrix M(0). According to ([2]), the covariance matrix M{t) of the time-evolved state p{t) 
satisfies 

^Mj^k{t) = i Tr (cjCki-i) [H, p]) = - Tr {p[H, CjCk]) for any j <k , (23) 

where we used the cyclicity of the trace. From this expression, it follows that the covariance 
matrix M{t) satisfies the equation 

|M(t) = [M(t),H] . (24) 

Eq. has the solution 

M{t) = R{t)M{0)R{tf where R{t) = e""* G S0{2N,R) . (25) 

The matrix R{t) can be computed in time 0{N'^) from the Williamson normal form of H. 
The latter can be computed in time 0{N^) (by diagonalizing H"^H), hence it follows that the 
evolution M i— )■ M{t) can be simulated in time 0{N^). Eq. f l25|) is sufficient for the purpose of 
simulating unitaries for fermionic linear optics computation starting from a Gaussian state. 

Let us discuss an alternative derivation of Eq. ( I25l) . which additionally provides a method 
for simulating the evolution of higher moments of a (possibly non-Gaussian) initial state p. It 
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also motivates the approach to dissipative dynamics discussed in Section 16.21 For this purpose, 
consider the adjoint action 



2N 



c,{t) = e^^*c,e-^^* = J2 mj,kCk (26) 

k=l 

of the unitary e~*^* on Majorana operators. Eq. fl2B]l can be shown as follows (cf. ^15^ Theo- 
rem 3]). Consider the derivative 

^ = [{-^)H,cM■ (27) 
Because [cjCk, q] = unless i G {j, k} and [cjCk, cj] = —2ck, we have 

[{-i)H{H), Cf\ = ^ iHm,„[cmC„, q] = - ^ iie,nCn ■ (28) 

n n 

Observe that if Cj{t) is a linear combination of the operators {q}, then so is ^Cj{t) because 
of fl27j) and fl28|) . Since this applies to t = 0, we conclude that Cj{t) is of the form specified 
on the Ihs. of fl2B]) for all t, that is, a linear combination of the operators {cfc} with some 
coefficients Rj^k{t)- It remains to find the matrix R{t). Rewriting (1271) in terms of R{t) and 
using ( |28ll . we obtain ^^^^ = — i?()f:)H by taking the anticommutator |{cfc, ■} on both sides, for 
k = 1, . . . , 2N. This shows that — H e so{2N, R) indeed generates R{t), proving the claim fl26p . 

Since the covariance matrix M{t) of the time-evolved state can be computed in the Heisen- 
berg picture as 

M,-fc(t) = Tr(e-'^Ve*^*c,Cfc) = Tr(pc,(t)cfc(t)) , 

the claim fl25l) immediately follows from fl26l) . For later reference, we point out that in this 
argument, we made use of fl26|) only to compute the product of two time-evolved operators, 
that is, in the form 

2N 

Cj{t)ck{t) = ^ R{t)j^iR{t)k,mCiCm ioT j ^ k . (29) 

It is also clear how this generalizes to higher moments: for example, if 

M,- fc,,(0) = Tr(p9CfeQ) (30) 

are the moments of a possibly non-Gaussian state p(0), the moments of the time-evolved 
state p(t) are given by 

Mj^k,i{t) = ^ Rjj'(t)Ri:^k'(t)Re/'(t)Mji^k',e' ■ 

j',k',l' 
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6 Simulating dissipative evolution 



In this section, we show how to simulate dissipative dynamics of the form ([3]). In more detail, 
since ([3]) preserves the Gaussian nature of a state p(0) according to Lemma [H it suffices to 
compute the covariance matrix M{t) of the time-evolved state p{t) = e^^{p). Our main result 
can be stated as follows: 

Theorem 1. Consider N fermionic modes and let C be a Liouvillian of the form ([3]). Let p 
he an even Gaussian state with covariance matrix M(0) G so(2A^, M). The covariance matrix 
M{t) of the time-evolved Gaussian state exp{tC){p) can be computed in time 0{N^). 

We give two different proofs of this statement: in Section I6.H we sketch how to solve the 
corresponding differential equation directly. This approach may be most numerically stable 
and relies on the Bartels-Stewart algorithm for finding the covariance matrix of a fixed point. 
Our second method, discussed in Section 16.21 is based on computing the Heisenb erg-evolved 
Majorana operators using the formalism of third quantization. This method can be adapted to 
simulate the evolution of higher moments starting from (possibly non-Gaussian) initial states. 

6.1 Simulation method based on evolution equation 

The definition of the covariance matrix M{t) implies that 

^Mj^k{t) = i Ti CjCkC{p{t)) = i Tr (a jCk)p{t) for any j < k. (31) 
The action of the adjoint Liouvillian £^ (cf. (|TOl) ) on observables can be written as 

£t(0) =^[H,0] + J2 LiP, L,] + [LJ , 0]L,. 

Choosing O ~ CjCk one can easily check that the commutators [H, O] and [O, L^] are quadratic 
and linear functions of the Majorana fermion operators respectively. It follows that preserves 
the subspace of operators spanned by the identity and CjCk with 1 < j < A; < 2N. Therefore 
Eq. (13T|) provides a closed linear differential equation that governs the time evolution of M{t). 
Simple algebra shows that 

—M(t) = XM(t) + M{t)yJ + Y, (32) 
dt 

where 

X = -H - 2(M + M*) and Y = 4z(M* - M) . 

This generalizes fl2^ . Recall that H is a real anti-symmetric matrix of size 2N x 2N parametriz- 
ing the quadratic Hamiltonian if, while M is a complex Hermitian matrix of size 2N x 2N 
parameterizing the Lindblad operators L^, see Section 14.11 It follows that X and Y are both 
real-valued 2N x 2A^-matrices, and Y = — Y-^ is antisymmetric. 
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Since any Lindblad equation has at least one steady statqj, there must exist at least one 
covariance matrix Mq which is a fixed point of Eq. f p2|) . that is, 

XMo + MqX^ + Y = 0. (33) 

Given any covariance matrix Mq satisfying Eq. fl33|) . the solution M{t) of Eq. f l32|) can be 
written as 

M{t) = Mo + e^*(M(0) - Mo)e^^*, t > 0. (34) 

Consider as an example a unitary evolution, that is, = for all /i. In this case M = 0, 
X = — H, and Y = 0. The fixed point covariance matrix can be chosen as Mq = and we 
recover fl25|) . 

Finding the fixed point covariance matrix Mq in the general case requires solving the system 
of linear equations Eq. fl55]) . where Mq is considered as an unknown vector of size N{2N — 1). 
A naive approach based on the Gaussian elimination would take time 0{N^) to find Mq. How- 
ever, this approach ignores the special structure of the problem. To the best of our knowledge, 
the most efficient method of solving the matrix equation Eq. (!33|) is the Bartels-Stewart algo- 
rithm [1] which has running time 0{N^). The key step of the algorithm is to perform a real 
Schur decomposition of X, that is, an orthogonal change of basis making X block upper trian- 
gular with blocks of size 1 and 2 on the main diagonal. In the new basis the resulting system 
of equations on matrix elements of Mq turns out to be quasi-triangular and can be solved in 
time 0{N^). For completeness, we sketch the Bartels-Stewart algorithm in Appendix Rl 

It is worth emphasizing that all matrices involved in Eq. f l5^ have bounded norm indepen- 
dent of t. This makes Eq. fl5^ suitable for numerical calculation of M[t). Indeed, taking into 
account that X + X^ = -4Re(M) < 0, one gets 

||gXi|| ^ ||g(X+xT)V2|| < 

Here the first inequality follows from Theorem IX. 3.1 of [3]. 

Finally, let us remark that Eq. f l32p can be solved directly even without knowing a fixed point 
Mq by transforming it into a homogeneous linear system. This method takes time only 0{N^), 
but unfortunately, it is computationally unstable for large evolution time t. Indeed, introduce 
an auxiliary 2A^ x 2A^ matrix K{t) whose time evolution is trivial, K[t) = K{0) = I. Then one 
can rewrite fl32|) as 

M{t) = XM(t) + M(t)X^ + Yis:(t), (35) 
k{t) = -X^K{t) + K{t)X.^, (36) 

with initial condition K{0) = I. Note that K(t) = / is indeed the only solution of Eq. (!36|) . 
Define a AN x 2A^ matrix Q such that 



M{t) 
Kit) 



^Note that for any integer n > there exists at least one (mixed) state p„ such that e^/"(p„) — p„. It 
imphes C{pn) = 0{l/n) for large n. Since pn belongs to a compact manifold, the sequence {p„} has a convergent 
subsequence. Its limiting point p is a (mixed) state that obeys C{p) = 0. 
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Then the system of Eqs. fl35|36p is equivalent to 

(lit) = Zn{t) + fi(t)X^, where Z = 

Its solution is 

n{t) = exp {Zt)n{0) exp (X'^t). (38) 

This can be computed in time 0{N^). Unfortunately, using Eq. (138!) for numerical calculations 
may be problematic due to the exponential growth of the factor exp (Zt). 

6.2 Simulation method based on 'third quantization' 

Without loss of generality, we assume that the even part = C(L) of the Liouvillian is 
specified by an antisymmetric matrix L G so(4A^, C) as in ([7]). We first derive the following 
analog of fl26l) (or, more precisely, fl29l) ) which applies to dissipative evolution. 

Lemma 2. For j < k, j, k E {1, . . . , 2N}, let A(j, k) be the AN x AN antisymmetric matrix 
with non-zero entries 

A(j, A;)2j_i+x,2fc-i+2; = for x,y E {0, 1} . (39) 

above the main diagonal and let A{k,j) = — A(j, /c). Furthermore, for p = 2j — l + x,q = 
2k — 1 + y , where x,y E {0, 1}, let T{p, q) be the antisymmetric 2N x 2N -matrix with non-zero 
entries 

T{2j -l+x,2k-l + y)j^k = {-if^^ for j < k and x,y E {0,1} (40) 
above the main diagonal, and let T{q,p) = —T{p,q). Then 

^ AN 2N 

exp{tCl){cjCk) = Y^J2Y1 WtfHj, k)R{t))r,sT{r, s)i^mCiCm for all j <k . (41) 

r,s=l l,m=l 

where R(t) = exp(L'l"t). 
Proof. We have by ([5]) and ([6]) 

^ 4N 

M = - 5^ A(j, k)p,,CpC,\I)) foTj^kE{l,..., 2N} (42) 

p,q=l 
^ AN 

CrCs\I)) = 2 5Z r^'^' for r 7^ s G {1, . . . ,4Ar} . (43) 

e,m=i 

Using (|9]) and unitality, we get 

exp{tCl)cpCg\I)) = Cp{t)c,{t) exp(t4)|/)) = Cp(t)c,(t)|/)) , (44) 



X 





Y 

-X^ 



(37) 
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where 

AN 

Cp{t) =exp{tCl)cpexp{-tCl) = ^R{t)p,rCr with R{t) = e^^^ e S0{4:N,C) (45) 

r=l 

is the Heisenberg time-evolved super- Major ana operator Cp. Equation (l45l) follows in the same 
way as (l26l) because we have 

[£(L),c,] = ^L,,fcCfc . (46) 

k 

in analogy with ( 128|) . Reinserting (H5!) into (jUj) gives 

exp(t£^)cpc,|/)) = J2 RpAt)RgAt)crCs\I)) ioi p ^ q . (47) 

Here we can restrict the sum to r 7^ s because R{t) is orthogonal and = I for all r. 

Equation p7|) is formally analogous to f l29|) . The claim follows by computing exp(t£^)|cjCfc)) 
using expression P2|) . linearity and P7|) . and then translating the result back with fH3|) . □ 

We can use Lemma [2] to compute the time-evolved covariance matrix M{t) of a Gaussian 
state p with covariance matrix M(0). Expressing the matrix elements using the Heisenberg 
picture gives 

M,- fc(t) = I TT{pexp{tCl){c,Ck)) for j <k (48) 

since exp(t£)(p) = exp(t£_|_)(p). Theorem [1] essentially follows by combining this equation with 
Lemma [21 However, showing that the resulting expressions can be evaluated in time 0{N^) 
requires some care. 

Let R = R{t) = exp(Ltt) be as in Lemma [2] and let A(j, k) and r(r, s) be defined by ( l39i) 
and (HOl) . respectively. With Lemma [2] and (HHl) we have 

^^•-'^W = E E ^)^)^'>«r(r, s),,„M,,„(0) = - Tr ((i?^A(j, A;)i?)fi(0)^) , (49) 

r,s=l £,m=l 

where fl{0) is the 4N x 4A^-matrix defined by the entries 

fi,,,(0) = Tr(r(r,s)M(0)^) . (50) 

Using the general identity Tr {{ABC)D'^) = Tr {{A'^DC^)B^) , we can rewrite as Mj k{t) = 
3^Tr((i^^](0)/^^)A(J,A;)^) or 

M,At) = ^Tr{n{t)A{j,kf) (51) 

where 

fi(t) = exp{'LH)n{0) exp{LHf. (52) 
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Clearly, Q{0) can be obtained from M(0) in time 0{N'^) since every entry ^2^,5(0) is a linear 
combination of a constant number of entries of M (cf. (1501) ). Similarly, the matrix M{t) can 
be computed in time 0{N'^) in an entrywise fashion according to fl3T]) given the matrix Q{t). 
Therefore, the claim of Theorem [T] follows since Q{t) can be computed from Q[0) in time 0{N^) 
using (152|) . To avoid exponentially growing terms arising from the expression exp(L^t) in fl52|) . 
we can decompose Q(t) as 

n{t) = n^it) + n^{t) (53) 

using the Jordan decomposition of L''', such that the matrix elements of Qb(t) are bounded for 
all t, whereas the matrix elements of Qu{t) are exponentially growing with t. Because M{t) is 
the covariance matrix of the time-evolved state exp(t£)(p), its entries are bounded for all t and 
we conclude from flSTl) that the contribution of flu{t) must vanish, i.e., 

M^,k{t) = ^TT{n,{t)A{j,kf) . 

To describe the decomposition ( I53ll in more detail, consider the Jordan normal form of L''', 

V = V^UK)V-' , (54) 

a 

where V is an invertible matrix. Here Ja{^a) is a Jordan block 

+ Nma 

where is the identity matrix of size nia x and Nm^ is the nilpotent matrix with ones 
on the first upper off-diagonal. With 

and exp(-J„(A„)t) = e"^"*^^^ (-t) we get 

= y ( 0e(''^-''')Xjt)V^-'f^(O)y5^,(-t) ) . (55) 

Clearly, all terms corresponding to pairs (a, (3) with Re{Xa — A^) > grow exponentially and 
must be absorbed in Qu{t)- Furthermore, for pairs with Re{Xa — A/3) = 0, the contribution of 
every nilpotent term grows polynomially with t. The remaining terms remain bounded, 
and we conclude that Qb{t) is given by 

\a,l3:Re{Xc-Xj3)<0 J 

+V [ e^^--'f^'>'lmaV-MO)VI^^ I . (56) 

\a,(S:Re{Xa-Xp)=0 J 
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Because computing the Jordan normal form fl5^ takes time 0{N^), the matrix Qb{t) (and 
hence also M(t)) can be computed in time 0{N^) using fl56l) . 

Finally, let us sketch the modifications required to simulate the evolution of higher moments 
such as (I5U1) starting from a non-Gaussian initial state. Clearly, Lemma [2] can be generalized to 
give an expansion of the Heisenberg evolved product exp(tC\^) (cjCkCi) in terms of a linear com- 
bination of Majorana monomials Cj'Ck'Cg/. Inserting this into Mj^k,e(t) = Tr(pexp(t£]^)(cjCfeQ)) 
immediately gives an explicit expression for the tensor of time-evolved moments M{t) in terms 
of the original moments M(0). 

A Computing the steady state 

In this appendix we sketch the Bartels-Stewart algorithm [1] for solving a matrix equation 

XM + MX^ + Y = (57) 

which determines the steady state covariance matrix M = Mq, see Eq. fl55]) in Section I^TTl Here 
X and Y are known real matrices such that Y"^ = — Y, while M is an unknown matrix to be 
found. All matrices have size nx n. As opposed to the third quantization method, where com- 
puting the steady state may involve manipulations with ill-conditioned invertible matrices, see 
Eq. the Bartels-Stewart method uses only transformations based on orthogonal (unitary) 
matrices. This is likely to offer a better numerical stability. 

The first step of the algorithm involves the Schur decompositiorH of X. This defines a 
unitary matrix U such that X = U'^IUJ is upper triangular. Introducing a new unknown 
matrix K = WMU we can rewrite Eq. fl57|) as 

XK + KX^ = Y (58) 

where Y = —WYU. Given any matrix Z, the j-th column of Z will be denoted Zj. Taking 
the n-th column of Eq. fl58p one arrives at 

(x + x;„/)ir„ = y„. (59) 

This is a triangular linear system of equations for the unknown vector Kn which can be solved 
in time O(n^). Taking the m-th column of Eq. (1581) one arrives at 

n 

{X + X*^^^I)Km = Yrn- Yl rn = n-l,...,2,l. (60) 

j=m+l 

Suppose the columns Km+i, ■ ■ ■ , Kn are already known. Then the right-hand side of Eq. (l60i) 
can be formed in time O(n^). We get a triangular system of equations for the unknown vector 
Km which can be solved in time O(n^). Proceeding inductively from m = n — 1 towards m = 1 
we can compute the entire matrix K in time O(n^). This gives a solution of Eq. fl57p . namely 
M = UKW. In general, this solution is neither real nor ant i- symmetric. However, one can 
easily check that matrices (M + M*)/2 and (M — M'^)/2 are solutions of Eq. fl57j) for any 
(complex) solution M. Hence we can always transform M into a real anti-symmetric solution. 

^To simplify notations, we slightly deviate from the original algorithm of [1] which adopted a real Schur 
decomposition. 
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Computing the Schur decomposition of a matrix X is a standard subroutine available in 
many numerical linear algebra tools (such as MATLAB). Theoretically, it can be computed in 
time 0{n^) by finding generalized eigenvectors of X and applying the Gram-Schmidt orthog- 
onalization. A more practical method used in pQ involves two steps. First, one transforms X 
to the upper Hessenberg form by a sequence of n — 2 Householder refiections, see [311 P- 346], 
which requires O(n^) elementary operations. Secondly, X is made upper triangular by the QR- 
algorithm. Each iteration of the QR-algorithm requires O(n^) elementary operations, however 
the required number of iterations is generally unknown. 
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